Data Load

library(dplyr)
library(knitr)
library(osmdata)
library(opentripplanner)
library(tidyverse)
library(sf)
library(ggthemes)
library(ggspatial)

Cary Public Parks data and NC state map load

CPP_parks <- st_read("https://data.townofcary.org/explore/dataset/parks-and-recreation-feature-map/download/?format=kml&timezone=America/New_York&lang=en")
opq(bbox = 'Cary NC USA') %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_xml(file = 'OTP/graphs/default/cary_streets.osm')
NC_state_plane <- "+proj=lcc +lat_1=34.33333333333334 +lat_2=36.16666666666666 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"

cary_street_features <- opq(bbox = 'Cary NC USA') %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_sf()

cary_streets <- cary_street_features$osm_lines %>%
  st_transform(crs = NC_state_plane)
ggplot(cary_streets) +
  geom_sf() +
  theme_map()

OTP data

path_otp <- otp_dl_jar("OTP")
path_data <- file.path(getwd(), "OTP")
path_otp <- paste(path_data, "otp.jar",sep = "/")

otp_build_graph(otp = path_otp, dir = path_data, memory =1024)
otp_setup(otp = path_otp, dir = path_data, memory =1024)
otpcon <- otp_connect()

Creating Isochrones

iso_10min_walk <- 
  otp_isochrone(otpcon = otpcon, fromPlace = CPP_parks, 
                mode = "WALK", cutoffSec = 600) %>%
  st_transform(crs = NC_state_plane) %>%
  mutate(mode = "walk")

iso_10min_drive <- 
  otp_isochrone(otpcon = otpcon, fromPlace = CPP_parks, 
                mode = "CAR", cutoffSec = 600) %>%
  st_transform(crs = NC_state_plane) %>%
  mutate(mode = "drive")

iso_all_modes <- rbind(iso_10min_drive, iso_10min_walk)

otp_stop()
right_side <- st_bbox(iso_all_modes)$xmax
left_side  <- st_bbox(iso_all_modes)$xmin
top_side <- st_bbox(iso_all_modes)$ymax
bottom_side <- st_bbox(iso_all_modes)$ymin

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = CPP_parks) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("By car", "By foot")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type = "cartolight", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = CPP_parks) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("By car", "By foot")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type = "stamenwatercolor", 
                      progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = CPP_parks) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_discrete(name = "Area that is reachable within 5 minutes",
                      labels = c("By car", "By foot"),
                      type = c("gray", "black")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type = "stamenbw", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = CPP_parks) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("By car", "By foot")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type = "hotstyle", 
                      progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = CPP_parks) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_discrete(name = "Area that is reachable within 5 minutes",
                      labels = c("By car", "By foot"),
                      type = c("gray", "black")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

ggplot(iso_all_modes) +
  geom_sf(data = cary_streets, color = "gray") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = CPP_parks) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("By car", "By foot")) +
  theme_map() 

Calcuate and compare isochrone areas

iso_areas <- iso_all_modes %>%
  mutate(area = st_area(iso_all_modes)) %>%
  st_set_geometry(NULL) %>%
  pivot_wider(names_from = mode, values_from = area) 

ggplot(iso_areas, 
       aes(x = as.numeric(walk), y = as.numeric(drive))) +
  geom_point() +
  scale_x_continuous(name = 
            "Area within a ten-minute walking distance\nof a public park\n(square km)",
            breaks = breaks <- seq(100000, 2000000, by = 100000),
            labels = breaks / 100000) +
  scale_y_continuous(name = 
            "Area within a ten-minute driving distance\nof a public park\n(square km)",
            breaks = breaks <- seq(0, 26000000, by = 2000000),
            labels = breaks / 100000) +
  theme_economist()